Abstract. We have investigated molecular distributions 
in protoplanetary disks, adopting a disk model with a 
temperature gradient in the vertical direction. The model 
produces sufficiently high abundances of gaseous CO and 
HCO + to account for line observations of T Tauri stars 
using a sticking probability of unity and without as- 
suming any non-thermal desorption. In regions of radius 
R > 10 AU, with which we are concerned, the tem- 
perature increases with increasing height from the mid- 
plane. In a warm intermediate layer, there are signifi- 
ed | cant amounts of gaseous molecules owing to thermal des- 
• orption and efficient shielding of ultraviolet radiation by 
the flared disk. The column densities of HCN, CN, CS, 
H2CO, HNC and HCO + obtained from our model are in 
\ good agreement with the observations of DM Tau, but 
are smaller than those of LkCal5. Molecular line profiles 
from our disk models are calculated using a 2-dimensional 
non-local-thermal-equilibrium (NLTE) molecular-line ra- 
diative transfer code for a direct comparison with obser- 
■ vations. Deuterated species are included in our chemical 
f^) ' model. The molecular D /H ratios in the model are in rea- 
sonable agreement with those observed in protoplanetary 

' disks. 
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1. Introduction 

It is well established from millimeter and infrared observa- 
tions that the birth of solar-mass stars is accompanied by 
the formation of a circumstcllar disk (Bcckwith & Sargent 
1996, Natta et al. 2000). These disks are important both 
as reservoirs of material to be accreted onto growing stars 
and as sites of planetary formation. Because the gas and 
dust in the disk are the basic components from which fu- 
ture solar systems are built, studies of their chemistry are 
essential to investigate the link between interstellar and 
planetary matter. Moreover, the chemical abundances and 
molecular excitation depend on physical parameters in the 
disks such as temperature and density, and on processes 
such as radial and vertical mixing. Thus, studies of the 
chemistry in protoplanetary disks can help to constrain 
their physical structure. 

Although CO millimeter lines arc routinely used to 
trace the gas and Keplerian velocity field in disks around 
classical T Tauri stars (e.g., Kawabe et al. 1993, Koerner 
et al. 1993, Dutrey et al. 1994, 1996, Koerner & Sargent 
1995, Saito et al. 1995, Guilloteau & Dutrey 1998, Thi 
et al. 2001), detections of other molecules are still rare. 
Dutrey et al. (1997) and Kastner et al. (1997) were the first 
to report observations of molecules such as HCO + , HCN, 
CN, HNC, H 2 CO, and C 2 H in the disks around GG Tau, 
DM Tau and TW Hya. Dutrey et al. found that the abun- 
dances of these species relative to hydrogen in the DM Tau 
and GG Tau disks arc lower than those in molecular clouds 
by factors of 5-100 (Dutrey et al. 1994, Dutrey et al. 1997, 
Guilloteau et al. 1999). The abundance ratio of CN/HCN, 
on the other hand, is significantly higher than in molec- 
ular clouds in all three disks. More recently, Qi (2000), 
van Zadelhoff et al. (2001) and Thi et al. (in preparation) 
have reported observations of molecules other than CO in 
the disks around the T Tauri stars LkCal5 and TW Hya, 
and in those around the Herbig Ae stars MWC 480 and 
HD 163296, confirming the trends of high CN/HCN ratio 
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and molecular depletion; the disk masses derived from CO 
(and its isotopes) are significantly smaller than those es- 
timated from the dust continuum assuming a gas to dust 
ratio of 100. 

Early models of the chemistry in disks considered 
mostly the one-dimensional radial structure in the cold 
midplanc (e.g. Aikawa et al. 1997, Willacy et al. 1998). 
Subsequently, it has been recognized that the vertical 
stratification of the molecules is equally relevant. Aikawa 
& Herbst (1999a) investigated the two-dimensional chem- 
ical structure within the so-called Kyoto disk model, rep- 
resentative of the minimum mass solar nebula. This model 
has low midplanc temperatures and is isothermal in the 
vertical direction. The molecular abundances were found 
to vary significantly with height Z from the midplane. 
At large radii (R >100 AU), the temperature is so low 
that most species, except H 2 and He, are frozen out onto 
the grains. This depletion is most effective in the mid- 
plane region (Z w 0) because of the higher density and 
hence shorter timescales for molecules to collide with and 
stick to the grains. In regions above and below the mid- 
plane, significant amounts of molecules can remain in the 
gas phase for longer periods because of the lower den- 
sities, and because of non-thermal desorption by cosmic 
rays and/or radiation (e.g., X-rays) from the interstel- 
lar radiation field and the central star. Aikawa & Herbst 
(1999a) suggested that the observed molecular line emis- 
sion comes mostly from this region. There is a height dis- 
tinction between stable molecules and radicals, however. 
In the surface region of a disk, radicals such as CN are very 
abundant because of photodissociation via ultraviolet ra- 
diation, whereas the abundances of the stable molecules 
such as HCN peak closer to the midplane. The molecular 
column densities obtained by integrating over height at 
each radius compare well with those derived from obser- 
vations of DM Tau (Dutrey et al. 1997), although detailed 
comparison through radiative transfer (i.e. calculation of 
molecular line intensities from model disks) is still lack- 
ing. The freeze-out of molecules in the midplane explains 
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the low average abundances of heavy-element-containing 
species relative to hydrogen, whereas the high abundance 
ratio of CN to HCN is caused by photodissociation in the 
surface layers. 

To obtain this good agreement with observations, 
Aikawa & Herbst (1999a) were forced to use the simpli- 
fying assumption that the probability S for sticking upon 
collision of a molecule with a grain is significantly smaller 
than unity. If the sticking probability is unity and if only 
thermal desorption is considered, the molecular column 
densities obtained in the Kyoto model are much smaller 
than observed, which suggests that either there is some 
efficient non-thermal desorption mechanism, or that the 
disk temperature is higher than assumed in the Kyoto 
model. Aikawa & Herbst (1999a) adopted an artificially 
low sticking probability S = 0.03 in order to reproduce the 
observed CO spectra, without specifying the non-thermal 
desorption process or modifying the temperature distribu- 
tion in the Kyoto model. Since adsorption is such a dom- 
inant process in the disk, the cause of the apparently low 
sticking probability should be considered more seriously. 

In order to explain the strong mid-infrared emission 
from disks, several groups suggested the possibility of 
higher dust temperatures than assumed in the Kyoto 
model due to efficient reception of stellar radiation by 
flared disks (e.g., Kenyon & Hartmann 1987, Chiang & 
Goldreich 1997, D'Alessio et al. 1998, 1999). In the two- 
layer Chiang & Goldreich (1997) model (C-G model here- 
after) , the upper layer (the so called "super-heated" layer) 
is directly heated by stellar radiation from the central 
star to temperatures T ;> 50 K at radii of ~ 100 AU. 
Also, recent observations of high-frequency lines (e.g. CO 
J = 6 — 5) support the possibility that the disk temper- 
ature is higher than assumed in the Kyoto model (Thi 
et al. 2001, van Zadelhoff et al. 2001). Willacy & Langer 
(2000) investigated the molecular distributions in the C-G 
model to see if the super-heated layer can maintain enough 
gaseous organic molecules to account for the observations. 
It was found, however, that molecules in this layer are de- 
stroyed by the harsh ultraviolet radiation from the star, 
whereas they are frozen out onto the grains in the cold 
lower layer. These authors therefore had to adopt a very 
high photodesorption rate in the lower layer to keep the 
molecules off the grains. 

In this paper we report an investigation of molecular 
distributions in another disk model with a vertical temper- 
ature gradient: the model of D'Alessio et al. (1998, 1999). 
These scientists obtained the temperature and density dis- 
tribution in steady accretion disks around T Tauri stars 
by solving the equations for local 1-D energy transfer (in- 
cluding radiation, convection and turbulence) and hydro- 
static equilibrium in the vertical direction. Whereas in the 
C-G model, the disk is divided into two discrete layers - 
super-heated and interior - the model of D'Alessio et al. 
gives continuous distributions of temperature and density. 
The differences in the temperature and density distribu- 



tions between these two models have a significant effect 
on the gaseous molecular abundances in the disk. Since a 
vertical distribution of density in the super-heated layer is 
not given explicitly in the C-G model, Willacy & Langer 
(2000) assumed a Gaussian density distribution with a 
scale height determined by the mid-plane temperature. In 
the model of D'Alessio et al., the gas is more extended 
than in a Gaussian distribution owing to the high tem- 
perature in the surface region. The higher densities (and 
thus the higher column densities) at large Z shield the 
lower layers from stellar ultraviolet radiation. In addition, 
the temperature variation in the vertical direction is more 
gradual in the model of D'Alessio et al. than the step 
function assumed in the C-G model. Therefore, compared 
with the C-G model, the model of D'Alessio et al. con- 
tains more gas in a warm and shielded layer, in which 
high abundances of gaseous molecules are expected. 

The rest of the paper is organized as follows. In §2 we 
describe the adopted model for protoplanetary disks and 
the chemical reaction network. Numerical results on the 
distributions of molecular abundances and column densi- 
ties are discussed in §3. In §4, molecular column densities 
and line intensities in the D'Alessio et al. model are com- 
pared with observations. Our conclusions and a discussion 
are given in §5. 

2. Model 

The disk model by D'Alessio et al. (1998, 1999) has been 
adopted in our work. In this model, the two-dimensional 
structure (i?, Z) of the disk is obtained by solving for 
hydrostatic equilibrium and energy transport in the in- 
direction. Various heating sources are considered, such as 
viscous dissipation of accretion energy, cosmic rays, and 
stellar radiation. Among them, stellar radiation is dom- 
inant at R > 2 AU, assuming typical parameters of T 
Tauri stars: T* = 4000 K, M* = 0.5 M , and R* = 2 R Q . 
We adopt a disk with an accretion rate M = 10~ 8 M 
yr^ 1 and viscosity parameter a = 0.01 as the fiducial, or 
standard, model, but also consider cases with the same 
a and differing accretion rates M = 10~ 7 M yr _1 and 
M = 10" 9 Mq yr" 1 (D'Alessio et al. 1999). The distri- 
butions of density and temperature in these three mod- 
els are shown in Fig. The surface density in the fidu- 
cial model is similar to that in the minimum-mass disk, 
which was adopted by Aikawa & Herbst (1999a), and is 
approximately an order of magnitude larger (smaller) in 
the model with the larger (smaller) M . In the fiducial disk, 
the masses inside 100 AU and 373 AU are about 0.017 M Q 
and 0.06 M Q , respectively. 

We have selected a few radial points from the model 
of D'Alessio et al. (viz. R = 26, 49, 105, 198, 290, and 373 
AU), divided the disk at each radius into several (30-40) 
layers, or slabs, depending on height Z 1 and calculated 
the molecular abundances in each layer as functions of 
time (Aikawa & Herbst 1999a). We have not included any 
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Fig. 1. The distributions of temperature T [K] (top row) and density n(H 2 ) [cm -3 ] (bottom row) in the D'Alessio et 
al. model with a viscosity parameter a = 0.01 are plotted as dotted lines. The top thick line in each panel shows the 
upper boundary of the model defined by P = 10~ 10 dyne cm~ 2 , which corresponds to a gas density of n(H 2 )= 1.4 10 4 
cm~ 3 with T — 50 K. From left to right the mass accretion rate is 10~ 7 , 10~ 8 , and 10~ 9 M Q yr -1 , respectively. In the 
top row, the grey scale shows the region in which the CO abundance relative to hydrogen is 10~ 6 — 10 (dark grey) 
and 10~ 5 — 10~ 4 (light grey). In the bottom row, the fractional HCO + abundance is 10 -11 — 10~ 10 (dark grey) and 
10- 10 - 10~ 9 (light grey). 



hydrodynamic motions in the disk, such as accretion or 
turbulence. The main goal of this paper is to investigate 
the effect of a vertical temperature gradient on molecu- 
lar abundances and line intensities. Although the model 
of D'Alessio et al. is an accretion disk model, the temper- 
ature distribution is determined by the irradiation from 
the central star and radiation transfer in the vertical di- 
rection, while the contribution of the accretion energy as 
a heat source is negligible in the region we are concerned 
with. In addition, we find that the chemical time scale is 
shorter than the accretion timescale, which is ~ 10 6 yr, 
in a large fraction of the molecular layers (section 3.2), so 
that molecular distributions obtained in this paper can be 
a reasonable approximation of reality. 

The chemical model and chemical reaction equations 
adopted in this paper are almost the same as those de- 
scribed in Aikawa & Herbst (2001). We use the "new 
standard model" for the gas-phase chemistry (Terzieva & 
Herbst 1998, Osamura et al. 1999), extended to include 
deuterium chemistry (Aikawa & Herbst 1999b). The ion- 
ization rate by cosmic rays is assumed to be the "stan- 
dard" value in molecular clouds, ( = 1.3 10~ 17 s _1 , be- 
cause the attenuation length for cosmic-ray ionization is 
much larger than the column densities in the outer regions 
of the disks, with which we are concerned (Umebayashi & 
Nakano 1981). Photoprocesses induced by ultraviolet ra- 
diation from the interstellar radiation field and from the 
central star are included. The ultraviolet flux from the 
central star varies with time and object, and can reach a 
value 10 4 times higher than the interstellar flux at R — 100 



AU (Herbig & Goodrich 1986, Imhoff & Appenzeller 1987, 
Montmcrle et al. 1993). This maximum value is adopted 
in this paper, as in Aikawa & Herbst (1999a). We assume 
that the ultraviolet radiation from the central star is not 
energetic enough to dissociate CO and H 2 . Self- and mu- 
tual shielding of H 2 and CO from interstellar UV is consid- 
ered as in Aikawa & Herbst (1999a). Chemical processes 
induced by X-rays from the central star are not included 
in this paper. 

Regarding gas-grain interactions, the surface forma- 
tion of H 2 , the surface recombination of ions and electrons, 
and the accretion and thermal desorption of ice mantles 
are included, but all other grain-surface reactions are not 
included. The sticking probability for accretion is assumed 
to be 1.0, unless stated otherwise. We have not considered 
any modifications to the grain-surface rate equation for 
H 2 formation (cf. Caselli et al. 1998), because, given our 
initial conditions, the rate of molecular hydrogen forma- 
tion is not important to the model. The total number of 
species and reactions included in our network are 773 and 
10446, respectively. 

The adopted elemental abundances are the so-called 
"low-metal" values (e.g., Lee et al. 1998, Aikawa et al. 
1999). The initial molecular abundances are obtained from 
a model of the precursor cloud with physical conditions 
n H = 2 10 4 cur 3 and T = 10 K at 3 10 5 yr, at which time 
observed abundances in pre-stellar cores such as TMC-1 
are reasonably well reproduced (Terzieva & Herbst 1998). 
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3. Results 

3.1. Vertical Distribution 

Figure 2 contains assorted vertical distributions at the 
outermost radius R = 373 AU in our fiducial disk model. 
The solid lines in Fig. |^ (a — c) show the physical param- 
eters: temperature, density, and A v . The attenuation of 
interstellar radiation (A* s ) is obtained from the equation 
N H 

Av = 1 r n -in 2 i =2 ma S> W 

1.59 lO^cm z 

where Nn is the vertical column density of hydrogen nu- 
clei from the disk surface to each point in the disk. The 
attenuation of stellar radiation (A^ tar ) is obtained via the 
same equation, but with ATh replaced by the column den- 
sity from the central star. Fig. |j (d — e) shows assorted 
molecular abundances at a disk age of t = 1.0 10 6 yr, 
which is typical for T Tauri stars. Significant amounts of 
molecules exist in the gas phase due to thermal desorp- 
tion at Z J> 100 AU, while most molecules are adsorbed 
onto grains below this height, where T < 20 K. As dis- 
cussed by van Zadelhoff et al. (2001), this molecular layer 
covers the region in which the lines of the main isotopes 
of the observed species become optically thick and thus 
where most of the observed emission arises. The results 
for models with different M are similar except that the 
height of the molecular layer is shifted in accordance with 
the distribution of A l J>, which is the main determinant of 
the vertical temperature distribution. 

The vertical distribution of the physical parameters 
(dashed lines) and molecular abundances in the Kyoto 
model are also shown in Fig. [2] (a — c, and /) for compari- 
son. The mass of the central star in the latter model is 0.5 
M©, as in the D'Alessio et al. model, so that the density 
distribution is modified from that adopted in Aikawa & 
Herbst (1999a). The sticking probability of neutral species 
onto grain surfaces is assumed to be 0.03, as in Aikawa & 
Herbst (1999a). It can be seen that the density distribu- 
tion in the D'Alessio et al. model is more extended than 
the Gaussian profile assumed in the Kyoto model (as well 
as in Willacy & Langer 2000), causing more efficient ultra- 
violet shielding of the warm layers just below the surface. 
Although we have not calculated the molecular distribu- 
tions in the C-G model, the width of its warm molecular 
layer at R = 373 AU can be estimated. In the C-G model, 
the Gaussian density distribution is similar to that in the 
Kyoto model, but the disk surface with A^ tar <j 1 mag 
is "super-heated" by stellar radiation. At R = 373 AU, 
the boundary between the super-heated layer and the in- 
terior region is located at Z ~ 280 AU, a value estimated 
from the Z— A^ tar relation in the Kyoto model (Fig. || 
c). Since the temperature in the interior region is lower 
than 20 K at this radius in the C-G model, warm gaseous 
molecules exist only at heights larger than 280 AU. The 
upper boundary of the molecular layer, which is deter- 
mined by photoprocesses, is estimated to be Z ~ 350 AU, 
again from the Kyoto model (Fig. || /). Therefore in the 



C-G model, the warm molecular layer, if any, is much nar- 
rower than in the D'Alessio et al. model. This estimate is 
consistent with the conclusion of Willacy & Langer (2000) 
that they need non-thermal desorption in the cold, more 
shielded layer to account for the observed molecular abun- 
dances within the C-G model. 

3.2. Radial Distribution of Column Densities 

Column densities are obtained by integrating the molec- 
ular abundances in the vertical direction. The column 
densities of assorted species as functions of disk radius 
are shown in Fig. || for three accretion disk models at 
t = 1.0 10 6 yr with different mass accretion rates. Stable 
neutrals such as H2CO and HCN show little dependence 
on radius, because these molecules are abundant only in 
regions with certain physical conditions and the mass con- 
tained in the layer with these physical conditions does 
not vary much with radius. For example, the HCN abun- 
dance is high (n(HCN)/n H ~ 10~ 10 - 10~ 9 ) only when 
n H £ 3 10 7 cm" 3 , T > 20 K, and A\ s £ 0.2 mag (see Fig. 
||) . The critical temperature of ~ 20 K is not the sublima- 
tion temperature of HCN, but that of CO, which is the 
dominant form of carbon in the gas phase. For HCN, at- 
tenuation of interstellar radiation is more important than 
that of stellar radiation because the interstellar radiation 
penetrates deeper into the disk due to the effect of geome- 
try. Radicals such as CN and C2H increase in column den- 
sity with radius because of the lower density and lower flux 
of the destructive stellar UV in the outer regions. Radical 
column densities are more sensitive than HCN to stellar 
UV, because their abundances peak at greater heights. 
Carbon monoxide is abundant (n(CO)/nn ~ 10 -4 ) in re- 
gions with T > 20 K and A l v s > 0.1 mag. The column 
densities of CO and HCO + abruptly change at R ~ 100 
AU, inside of which the temperature in the midplane is 
higher than 20 K. These characteristics of the radial dis- 
tribution arc similar to those in the Kyoto model (Aikawa 
et al. 1996, Aikawa & Herbst 1999a). 

The amount of gas existing under physical conditions 
conducive for large abundances of molecules does not vary 
significantly among the three disk models with different 
accretion rates (and thus with different disk mass), either. 
Therefore, most (but not all) molecular column densities 
vary only by a very small factor among the three disk 
models, even though the total (H2) column density varies 
by two orders of magnitude. In the model with a mass ac- 
cretion rate of 1.0 10 -9 M Q yr -1 , the region with density 
nH ~ 10 5 — 10 6 cm -3 , at which CN is abundant, is more 
shielded from stellar UV (Fig. [I]), and thus the CN column 
density is higher than in the other two models. 

Molecular column densities in our fiducial disk model 
at an earlier time of t = 1.0 10 5 yr are shown with thick 
dot-dashed lines in Fig. ^. The variation in column den- 
sity for most molecules during 10 5 — 10 6 yr is less than 
a factor of 2; two exceptions are CS and H2O. In regions 
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Fig. 2. Vertical distributions at R — 373 AU of (a) temperature, (b) density (rtjj = 2n(H2) + n(H)), (c) attenuation 
of the interstellar radiation (A* s ) and stellar radiation (A^ r ), (d — e) molecular abundances in the D'Alessio et al. 
model, and (/) molecular abundances in the Kyoto model with S = 0.03. In panels (a — c), the physical parameters 
of the D'Alessio et al. and Kyoto models are shown via solid and dashed lines, respectively. In panels (d — /), the disk 
age is assumed to be t = 1.0 10 6 yr. 
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Fig. 3. The column densities of assorted molecules as functions of radius in models with mass accretion rate 1.0 10 
(solid lines with crosses), 1.0 10~ 8 (solid lines with closed circles), and 1.0 10 -9 (dashed lines with open circles) M 
yr" 1 . The viscosity parameter is fixed at a = 0.01, while the disk age is t = 1.0 10 6 yr. The column densities at 
t = 1.0 10 5 yr in the model with mass accretion rate 1.0 10 yr -1 are shown with thick dot-dashed lines. 



with A l J> smaller than a few mag, which covers a large 
fraction of the molecular layer, the chemical timescales 
are short (<; 10 5 yr) because of photoprocesses. In the 
more shielded portion of the molecular layer at smaller Z, 
S-bearing gaseous molecules decrease in abundance after 
10 5 yr, since most sulfur is adsorbed onto grains in the 



form of CS, SO and OCS. Similarly, H2O gas decreases at 
~ 10 6 yr, because most oxygen which is not in CO is ad- 
sorbed as H2O ice. On the other hand, abundances of other 
C-bearing molecules reach pseudo-steady-state values in a 
relatively short timescale (<; 10 4 yr), and do not show sig- 
nificant time variation during 10 — 10 6 yr in the more 
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shielded region because CO gas, their chemical precursor, 
remains the dominant component of carbon for more than 
1 10 6 yr. The pseudo-steady-state gas-phase chemistry of 
non-volatile carbon-containing species is balanced for a 
considerable period by formation reactions starting from 
CO and depletion onto the dust particles. 

Values of the sticking probability S are estimated to lie 
in the range 0.1 <J S <; 1.0 (Williams 1993, and references 
therein). In order to check the dependence of the molecu- 
lar column densities on S, we have performed calculations 
with S = 0.1 in addition to our fiducial value of unity. 
Column densities of radicals such as CN and C2H do not 
depend on S, because they are abundant in the surface 
layer, in which adsorption is not the dominant process. 
Among the more stable species, some show significant de- 
pendence on S; the column densities of CO2 and OCS are 
larger by an order of magnitude in the model with the 
lower sticking probability at R = 373 AU and t = 1 10 6 
yr. But the effect of S is smaller for many other species; 
at R = 373 AU and t = 1 10 6 : for example, the column 
densities of CS, HCN, and H2CO are larger only by fac- 
tors of 3.2, 1.7, and 1.3, respectively, in the case of lower 
5*. There are two reasons for this small dependence on 
S. First, adsorption is not always the dominant process in 
the molecular layer, depending on species and height from 
the midplane. Second, the higher abundance of gaseous O2 
in the case of lower S reduces the abundance of the car- 
bon atom, and thus reduces the formation rate of organic 
molecules, which counteracts the lower adsorption rate. 

3.3. D/H Ratios in Molecules 

In the disks around LkCal5 and TW Hya (Qi 2000, Thi 
et al. in preparation), the deuterated species DCN and 
DCO + have been detected, and the ratios of deuterated 
to normal species (DCN/HCN and DCO+/HCO+) are es- 
timated to be ~ 0.01. The deuterated species HDO has 
been detected in LkCal5, but there is no estimated ratio 
of HDO to normal water. Although the estimated ratios 
might include uncertainties as discussed in §4, they would 
not alter the important conclusion that the molecular D /H 
ratios are higher than the cosmic elemental abundance ra- 
tio (D/H« 1.5 10~ 5 ) by orders of magnitude. It is inter- 
esting to see if our model reproduces the high D/H ratios 
and if these ratios can be used as probes of the chemical 
processes and physical conditions in disks. Fig. |] shows 
some column density ratios of deuterated species to nor- 
mal species in our models as functions of radius. The ratios 
are similar to the observed values of 0.01 for both HCN 
and HCO + . The mechanism of deuterium enrichment in 
disks is similar to that in molecular clouds; due to en- 
ergy differences between deuterated and normal species, 
molecules such as H^" have a high D/H ratio, which prop- 
agates to other species through chemical reactions (Millar, 
Bennet & Herbst 1989, Aikawa & Herbst 1999b, Roberts & 
Millar 2000). The D/H ratios decrease as the temperature 
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Table 1. Calculated (t = 1.0 10 6 yr) molecular column 
densities (cm -2 ) at R — 373 AU compared with observa- 
tions 



Species M [MQyr -1 ] DMTau a LkCal5 

10~ 7 10~ 8 interferometer 13 single dish c 

H 2 1.5(24) d 1.3(23) 4.3(21) 3.8(21) 

CO 1.1(18) 1.1(18) 5.8(17) 5.7(16)° 1.6(18)° 9.0(17)* 

HCN 2.1(12) 2.3(12) 1.5(12) 2.1(12) 0.02 - 1.2(15) g 7.8(13) 

CN 3.8(12) 4.7(12) 1.7(13) 9.5-12(12) 9.1-25(13) 6.3(14) 

CS 4.9(11) 5.6(11) 1.4(12) 6.5-13(11) 1.9-2.1(13) 2.2(14) 

H 2 CO 2.9(12) 3.3(12) 3.8(12) 7.6 - 19(11) 3.0 - 22(13) 

HCO+ 9.0(12) 9.9(12) 4.8(12) 4.6 - 28(11) 1.5(13) 1.4(13) 

C 2 H 6.2(12) 6.0(12) 7.9(12) 4.2(13) 

HNC 2.0(12) 2.3(12) 1.5(12) 9.1(11) < 5.4(12) 

OCS 3.1(10) 2.8(10) 5.0(9) < 2.9(13) 

CH 3 OH 6.4(8) 7.1(8) 6.6(8) 7.3 - 18(14) < 9.4(14) 

DCN 5.5(10) 6.4(10) 3.7(10) 1.0(13) 

HDO 2.1(12) 2.4(12) 5.7(11) 2.3 - 6.8(14) 

N 2 H+ 1.9(12) 1.9(12) 6.0(9) < 7.6(11) < 5.7(12) < 5.9(13) 



a Derived from single-dish data by Dutrcy et al. (1997) (sec text). 

b Derived from interferometer data by Qi (2000). The values in this column do 

not necessarily refer to 373 AU (sec text). 
° Derived from single-dish data by Thi et al. (in preparation) assuming a disk 

radius of 100 AU (sec text). 
d 0(6) means a 10 b . 

° Estimated from C 1S assuming C 16 0/C 1S 0= 500. 
f Estimated from 13 CO assuming 12 CO/ 13 CO= 60. 

g Lower value is estimated from H 12 CN and higher value is an upper limit 
estimated from H 13 CN assuming H 12 CN/H 13 CN= 60. 



rises, because the energy differences are less significant at 
higher temperatures. Thus, the D/H ratio decreases in- 
wards (Fig. ||), because the temperature in the gaseous 
molecular layer is higher at inner radii. Because the tem- 
perature in the molecular layer does not depend much on 
the disk mass (see Fig. |l|), the differences between models 
with different M are small. 

The ratio of HDO/H2O shows a more complicated 
behavior than described above. It has a local peak at 
R ~ 100 AU, and the peak is higher in the more mas- 
sive disk model. At R ~ 100 AU, the midplane is at 14-16 
K for the three disk models, at which temperature CO is 
almost completely frozen onto grains but O can marginally 
be kept in the gas phase to produce water vapor. The D/H 
ratio in the vapor is enhanced by the CO depletion in the 
region, since CO is one of the main destruction channels of 
H2D+ (Brown ct al. 1989). The layer at this critical tem- 
perature (14 K<; T <J 16 K) is thicker in the more massive 
disk. At smaller radii, the midplane temperature is higher, 
which lowers the D/H ratio. At larger radii for disks with 
M = 10~ 7 and 10~ 8 M Q yr -1 , the midplane temperature 
is so low that O cannot remain in the gas phase to pro- 
duce water vapor. The abundance of HDO has a sharp 
peak in a thin layer of 14 K <J T <J 16 K offset from the 
midplane (Fig. || d). In the outer disk of the M = 10 -9 
Mq yr -1 model, the temperature is slightly higher than 
20 K even in the midplane (see Fig. [I]), which causes a 
lower HDO/H2O ratio than for the other two models. 

4. Comparison with observation 
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Fig. 4. Column density ratios of deuterated species to normal species as functions of radius in the three accretion disk 
models at an age of t = 1.0 10 6 yr. The models are represented by the same lines as in the previous figure. 
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Fig. 5. Molecular line profiles from LkCal5 (solid line), and the D'Alessio et al. model with accretion rates of 1.0 10 -9 
(dot-dashed line), 1.0 10~ 8 (thick solid line), and 1.0 10 -7 (dashed line) M Q yr -1 . The dotted line in the HCN profile 
represents the disk with accretion rate of 1.0 10~ 8 M© yr" 1 , but HCN abundance at each position in the disk is 
artificially enhanced by a factor of 10 from our calculated values. Similarly, in the CN profile abundance is enhanced 
by a factor of 50 for the dotted line. The disk radius, age, and inclination are 400 AU, 1 10 6 yr and 60 degrees, 
respectively. 



4-1- Column densities 

In the subsequent paragraphs, we discuss comparisons of 
our calculated column densities for assorted species at a 
particular radius with both single dish and interferometric 
data. It must be recognized, however, that it is not really 
possible to derive reliable columm densities at a particular 
radius from unresolved single dish data, making much of 
the subsequent discussion less quantitative than would be 
desirable. 

Table compares calculated molecular column densi- 
ties obtained with three different accretion rates at a time 
of t = 1.0 10 6 yr and a radius of R — 373 AU with those 
estimated from the observations of DM Tau (Dutrey et 
al. 1997) and LkCal5 (Qi 2000, Thi et al. in preparation). 
In general, it is difficult to estimate the total H 2 column 
density in the disk, especially at the outer radius. Dust 



observations suffer from uncertainties in the dust opac- 
ity at millimeter wavelengths, which depends on the grain 
size distribution. Moreover, the dust continuum is very 
weak at the outer radii and difficult to detect. The H2 
column density cannot be directly estimated from molec- 
ular lines, because the molecular abundances relative to 
hydrogen are not known. Therefore Dutrey et al. (1997) 
estimated the H2 column density and averaged molecular 
abundances for DM Tau in a different way, paying atten- 
tion to the critical density for excitation of the molecular 
lines and their optical depth. For the DM Tau disk, inter- 
ferometric observations of 12 CO (J = 2 — 1) show that the 
CO gas extends to ~ 800 AU from the central star. Com- 
bining the different constraints, they obtained an H2 den- 
sity n(R,Z) = 5 10 5 (i?/500AC7)- 3 exp[-(Z/i?) 2 ] cm" 3 
over this range, in which H is the scale height of the disk, 
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with H « 175 AU at R = 500 AU. From the line intensi- 
ties obtained in single dish telescopes, they subsequently 
derived molecular abundances with respect to hydrogen 
assuming that the abundances are constant over the entire 
disk. We obtained the molecular column densities of DM 
Tau in Table |l| by vertically integrating this disk model 
using the average abundances listed in Table 1 of Dutrey 
et al. (1997). 

Although the DM Tau disk extends to about 800 AU 
from the central star, and the outer radius would have a 
larger contribution to the line intensities because of the 
larger surface area, we list the calculated column density 
at R = 373 AU, since the D'Alessio et al. model extends 
only out to this radius (the H 2 column density at R = 800 
AU is about 3 times less than that at R = 373 AU in the 
DM Tau disk). The model with an accretion rate M = 
10~ 9 M Q yr~ x reproduces reasonably well the vertical H2 
column density of DM Tau, and agrees with other DM 
Tau observations to within a factor of 2, except for CO 
and C2H. The former is overestimated, while the latter is 
underestimated. It should be noted that we have assumed 
that the stellar UV is not hard enough to dissociate CO, so 
that our model may overestimate the CO column density. 
Inclusion of CO photodissociation by stellar UV might 
also increase C2H and CN abundances, since more carbon 
is released from CO in the photodissociation region, in 
which C2H and CN are mainly formed (van Zadelhoff et 
al. in preparation). 

Qi (2000) performed interferometric observations on 
LkCal5 and estimated beam averaged molecular column 
densities based on the velocity-integrated intensity mea- 
sured over a much smaller beam than used to determine 
the DM Tau abundances. The vertical column densities 
are lower than the listed values by a small factor, which 
is less than 2 if the inclination is <; 60 degrees. Since the 
beam size is about 0.6" — 13" depending on the frequency 
of the line, the estimated values do not necessarily refer to 
373 AU. If the emission is resolved and if we assume that 
the molecular column densities do not vary much with ra- 
dius (as indicated by our disk models for R > 100 AU), 
we can take the values listed by Qi (2000) to apply to 373 
AU. The typical beam size of Qi (2000) is 3" - 4" (~300 
AU radius at the distance of LkCal5), so the majority of 
the results will be just resolved, making it reasonable if 
not perfect to compare the observations with our result at 
R = 373 AU. 

In addition to the interferometric work on LkCal5, Thi 
et al. (in preparation) derived molecular column densities 
for this disk from high-frequency single dish observations 
of various species. Their beam-averaged values obtained 
with their assumption that the disk has a radius of 100 AU 
are included in Table 1 and differ from the interferometric 
column densities by up to an order of magnitude. Both 
values are much higher than the column densities found 
for DM Tau, so the agreement between our model and 
observations is worse for LkCal5 than for DM Tau. For 



CO and HCO + the difference is less than a factor of 2, but 
the column densities of other species are about 1-2 orders 
of magnitude higher than the model results for all three 
mass accretion rates. Methanol in LkCal5 appears to be 
a special case; its calculated column density is more than 
five orders of magnitude too low, presumably because it is 
produced by the hydrogenation of CO on grain surfaces, 
which is not included in our model, and then evaporated 
into the gas. 

The use of 100 AU for LkCal5 as a reference point by 
Thi et al. is somewhat arbitrary. If the values obtained 
by Thi et al. are taken to refer to a 373 AU radius disk, 
their LkCal5 column densities need to be reduced by a 
factor (100/373) 2 , and become closer to those found for 
DM Tau. Thi et al. and Qi (2000) also present data for 
a few other disks (TW Hya, HD 163296, and MWC 480) 
and determine column densities and abundances in a con- 
sistent way. Indeed, the DM Tau column densities and 
abundances are generally lower than the values for other 
disks, whereas those for LkCal5 are among the highest. 
Thus, these two sources appear to bracket the range of 
observed values. 

4-2. Line intensities and profiles 

In addition to the comparison with estimated molecular 
column densities at a single radius, a more direct com- 
parison via line intensities from the entire model disk has 
been made. Since the most complete single-dish and inter- 
ferometer data set is available for LkCal5, we restrict our 
efforts to this source. In estimating the molecular column 
densities from the observations, Qi (2000) assumed local 
thermal equilibrium (LTE) with an excitation tempera- 
ture of 40 K throughout the LkCal5 disk, thereby deriving 
a mean column over the beam. In this paper, however, we 
have shown that a temperature gradient is important for 
the characteristics of the gaseous molecular layer in the 
disk and that the abundances are strongly varying with 
R and Z. Also, the excitation of the molecules and the 
line emission depend on the density structure; van Zadel- 
hoff et al. (2001) have shown that the assumption of LTE 
is not always valid for high frequency lines of molecules 
with high critical densities. Hence, it is better to compare 
directly simulated line emission with observations. 

We calculate here the excitation of the molecules us- 
ing the 2-dimensional (2D) NLTE molecular line radia- 
tive transfer code of Hogerheijde & van der Tak (2000). 
From the resulting level populations, the line profiles can 
be computed, taking into account the inclination of the 
source. NLTE molecular line radiative transfer in more 
than one dimension has been used in star-formation re- 
search only recently (e.g., Park & Hong 1995, Juvela 
1997). The need for a full treatment of the radiative trans- 
fer in disks follows from the non-locality of the problem. 
The level populations depend both on the local param- 
eters (temperature, density and radiation field) and on 
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the global radiation field, which in turn depends heavily 
on the optical depth of the medium. The main problem is 
the slow convergence of the level populations. The adopted 
code by Hogerheijde & van der Tak (2000) uses an accel- 
erated Monte Carlo method to enhance convergence. A 
more elaborate discussion of the methods is given in van 
Zadelhoff et al. (2001). 

The main input parameters to the radiative transfer 
calculations in addition to the molecular data (collision 
rates, Einstein A coefficients) are the 2D distributions of 
temperature, density and abundance of the molecule of 
interest at each position (R, Z) . The latter is taken di- 
rectly from the chemical models at a chosen age. Other 
important parameters are the turbulent width, taken to 
be 0.2 km s _1 (van Zadelhoff et al. 2001), the systematic 
Keplerian velocity field for an assumed stellar mass, the 
thermal line width, the size of the disk, and the inclina- 
tion of the source. Qi (2000) estimates a disk inclination 
of ~ 58 ±10 degrees, and an outer radius of the CO disk of 
435 AU for LkCal5. In the following, we assume a disk in- 
clination of 60 degrees and a size (outer radius) of 400 AU. 
The dependence of line emission and profiles to disk radius 
and inclination can be roughly estimated (Omodaka et al. 
1992, Beckwith & Sargent 1993), e.g., for optically thick 
species the peak intensities of molecular lines are propor- 
tional to the square of the disk radius. Dependence of the 
integrated intensities on disk inclination is also discussed 
in van Zadelhoff et al. (2001), and is found to be less than 
a factor of two when the inclination is varied from to 60 
degree. 

Fig. H compares simulated line profiles for CO (J = 
6 - 5, 3 - 2, and 2 - 1), HCO+ (4 - 3), HCN (4 - 3), 
and CN (3 — 2) from the fiducial D'Alessio et al. model 
(thick solid line) with lines from LkCal5 observed with 
single-dish telescopes by van Zadelhoff et al. (2001) (solid 
line). Line profiles based on the D'Alessio et al. model 
with other accretion rates - M = 10~ 9 (dot-dashed line) 
and M = 10~ 7 (dashed line) M yr 1 - are also shown. 
In spite of the fact that the column densities of CO and 
HCO + in our model are slightly smaller than those es- 
timated by Qi (2000), the calculated intensities of these 
species are higher than observed in LkCal5 by a factor of 
2 — 3, which is caused by the higher disk temperatures in 
our model. The vertical temperature gradient lowers the 
optical depth of the disk, which further enhances emis- 
sion intensities (van Zadelhoff et al. 2001). As opposed to 
CO and HCO+, the model intensities of CN and HCN are 
much lower than those observed in LkCal5. Those lines 
are optically thin in our model, and about 10 times more 
HCN and at least 50 times more CN are needed to fit the 
observed profiles, which is consistent with the compari- 
son of column densities. The dotted lines in the CN and 
HCN panels show model profiles in which the molecular 
abundances are artificially enhanced. We conclude that 
CN and HCN in LkCal5 are much more abundant than 
in our model. 



4-3. Discussion 

What are the causes of the disagreement between our 
model results and LkCal5, and the difference between DM 
Tau and LkCal5 ? Let us first discuss the uncertainties in 
the estimate of the molecular column densities and size of 
the emission region. Although there are some differences 
in observed intensities, the lines in LkCal5 are, in fact, not 
much stronger than those in DM Tau, at most a factor of 
2-3 after scaling the data to the same beam. The main 
reason for the different column densities in DM Tau and 
LkCal5 then seems to be the adopted size of the emission 
region, because the column densities are inversely propor- 
tional to the square of the size (radius) of the emission 
region. For CO, where interferometric observations have 
been performed both for DM Tau and LkCal5, the DM 
Tau disk is found to be somewhat larger than the LkCal5 
disk — DM Tau is about 800 AU in radius, while LkCal5 
is 435 AU. But whether this larger emission region in 
DM Tau holds for other species is not yet known without 
some degree of ambiguity because the DM Tau data come 
mainly from single-dish observations, and the size of the 
region and average molecular abundances are estimated 
from a model with constant abundances throughout the 
disk. For LkCal5, in contrast, the sizes of the emission 
regions are directly derived from interferometric observa- 
tions and are at most a few arcseconds. Assuming that the 
interferometer resolves the emission, these column densi- 
ties should therefore not be affected by uncertainties in 
the size of the emission region. Although interferometry 
is more direct, and thus seems more reliable for obtaining 
the size of the emission region, it should be noted that the 
derived size depends on the signal to noise ratio and dy- 
namic range achieved in the observations. In fact, the size 
of the emission region obtained via interferometry seems 
to contain uncertainties; the size of the CO ( J = 2 — 1) 
emission region around LkCal5 is estimated to be ~ 600 
AU by Duvert et al. (2000), which is larger than the 435 
AU obtained by Qi (2000). 

Since our 2D modeling procedure convolves the cal- 
culated molecular emission with the actual beam of the 
observations, a larger disk size cannot explain the discrep- 
ancy with the LkCal5 interferometer observations, but it 
does affect the calculated single dish intensities. Such an 
increase would not be sufficient to explain the discrepan- 
cies for LkCal5, however. For example, if we assume an 
outer disk radius of 600 AU, the calculated HCN and CN 
line intensities from our model disk (Fig. |J) are increased 
by a factor of 2 at most, while the calculated CO and 
HCO + lines become even stronger compared with the ob- 
served profiles. There are a few possible explanations for 
these discrepancies. First, our model might overestimate 
the CO column density and underestimate the radical col- 
umn density because we do not consider dissociation of CO 
via stellar UV radiation, as mentioned above. Detailed 
consideration of stellar UV radiation, including the CO 
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and H2 dissociation, will be reported in a forthcoming pa- 
per. Indeed, CN is enhanced by an order of magnitude de- 
pending on the treatment of the radiation field, although 
HCN is not much changed. The inclusion of X-rays might 
be another solution. X-rays cause ionization and dissocia- 
tion, which enhance chemical activity, and hence increase 
the transformation of CO to other organic molecules such 
as CN and HCN (Aikawa & Herbst 1999a, 2001). Different 
X-ray fluxes might also account for the differing molecu- 
lar column densities in DM Tau and LkCal5, if they are 
intrinsic. X-rays also cause non-thermal desorption, which 
might enhance the CN and HCN abundances in the gas 
phase (Najita, Bergin & Ullom 2001). Finally, of the disks 
around T Tauri stars surveyed so far, LkCal5 stands out 
as the disk with the strongest molecular lines and richest 
chemistry in the interferometer and single-dish data (see 
§4.1, Qi 2000, Thi et al. in preparation). 

5. Conclusion and Discussion 

We have investigated the molecular distributions in proto- 
planetary disks by combining the "new standard" chem- 
ical model with the physical model of D'Alcssio et al., 
which has a temperature gradient in the vertical direction. 
The calculated molecular column densities are in reason- 
able agreement with those estimated from single-dish ob- 
servations of the DM Tau disk without the assumptions in 
previous calculations of non-thermal desorption and/or an 
artificially low sticking probability. In the warmer inter- 
mediate layers of our disk models, there are large amounts 
of gaseous molecules owing to thermal desorption and ef- 
ficient UV shielding caused by large gas densities at large 
heights from the midplane, a phenomenon known as flar- 
ing. Gaseous molecules are abundant only in regions with 
certain physical conditions. The volume of the layers with 
these conditions, and thus the column densities of gaseous 
molecules, are not proportional to the total (H2) column 
density. Column densities of abundant molecules such as 
CN, HCN and HCO + do not vary by more than a factor 
of three during the period t <~ 10 5 — 10 6 yr. Sulfur-bearing 
molecules and H 2 show larger temporal variations. 

Comparison of our model results with those of Willacy 
& Langer (2000), who adopted the C-G disk model with 
a Gaussian density distribution, indicates that gaseous 
molecular abundances are sensitive to the vertical struc- 
ture of the disk model; in their model, molecules in the 
super-heated upper layer are destroyed by the harsh ultra- 
violet radiation from the star, while sufficient UV shield- 
ing is available in the warm upper layers of the D'Alessio 
et al. model. Deuterated species are also included in our 
chemical model. The molecular D/H ratios we obtain are 
in reasonable agreement with those observed in protoplan- 
etary disks. 

Despite our agreement with observations of DM Tau, 
the molecular column densities obtained in our models 
are smaller than those observed around LkCal5, except 



for CO and HCO + . The estimated column densities of all 
observed molecules around LkCal5 are higher than those 
around DM Tau by roughly an order of magnitude. This 
difference in derived molecular column densities in the two 
objects seems to derive, at least partially, from different 
and/or uncertain sizes of the emission regions. Compari- 
son with other sources shows that some of the difference 
is likely to be intrinsic, and other physical parameters or 
processes, such as X-ray ionization and dissociation, are 
needed to account for the high column densities in the 
LkCal5 disk. 

In addition to the calculation of molecular abundances 
and column densities, we have solved the equation of radi- 
ation transfer to obtain line profiles from our model disks, 
which can be directly compared with observations. Such a 
comparison is a much more detailed test of theory than is 
a comparison of column densities, since line intensities de- 
pend not only on the molecular column densities but also 
on the density and temperature of the molecular layer and 
the variation of the abundance of the molecule with R and 
Z. The line intensities of HCN and CN obtained from the 
theoretical models are lower than the observed intensities 
in LkCal5, as expected from the comparison of column 
densities. 

There are still several uncertainties in the vertical 
structure of protoplanetary disks, which might affect our 
results. Firstly, Chiang & Goldreich (1997) argue that the 
gas temperature could be lower than the dust temperature 
in the upper layers; although an important heat source of 
the gas is collisions with super-heated dust particles, gas- 
dust collisions are not frequent enough to equilibrate the 
gas temperature and grain temperature because of the low 
density. With a lower gas temperature, the disk would be 
less flared than in the model of D'Alessio et al., which 
assumes equal temperatures of gas and dust. Glassgold & 
Najita (2001) have pointed out, however, that gases in the 
upper layer can be heated by X-rays, which were not con- 
sidered by Chiang & Goldreich (1997). Since Glassgold & 
Najita (2001) have listed the temperature only in the in- 
ner radius (R = 1 AU), we made a rough estimate of the 
disk surface temperature for the outer radii R ~ 100 — 300 
AU based on the work of Maloney et al. (1996), which sug- 
gests that the surface temperature of the X-ray irradiated 
disk is indeed higher than the midplane (interior) temper- 
ature of the C-G and Kyoto models. Moreover, the UV 
photons can heat the gas through the photoelectric effect 
on grains and PAHs, as in models of photon-dominated re- 
gions. Hence we can at least expect disks to fall off more 
slowly in density with increasing height than in a simple 
Gaussian distribution, although detailed studies on the 
heating and cooling balance between gas and dust are de- 
sirable in order to obtain an accurate vertical structure for 
the disk. Another uncertainty lies in the size and distribu- 
tion of dust particles. D'Alessio et al. (2001) and Chiang 
et al. (2001) have noted that their original models are ge- 
ometrically too thick compared with the observations of 
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cdgc-on disks, which suggests dust sedimentation and/or 
growth in the disk. Because the molecular abundances in 
our model depend on the efficiency of UV shielding by 
"small" (i.e. interstellar) dust grains (Aikawa & Herbst 
1999a), we might have overestimated these abundances. 
Although a more detailed approach with dust sedimenta- 
tion is beyond the scope of this paper, we emphasize that 
molecular abundances can help to resolve uncertainties in 
dust evolution and disk structure. 
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